EPJ manuscript No. 

(will be inserted by the editor) 



Force distributions and force chains in random stiff fiber networks 



' Claus Heussinger 1 and Erwin Frey 1 

| Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience (CeNS), Department of Physics, Ludwig- 
s ! ■ Maximilians-Universitat Miinchen, Theresienstrasse 37, D-80333 Miinchen, Germany 



o 
O 

(N 



<4-l 

o 

c3 



T3 
C 

o 
o 

(N 
> 

in 

in 
o 

o 



X 



February 1, 2008 

Abstract. We study the elasticity of random stiff fiber networks. The elastic response of the fibers is 
characterized by a central force stretching stiffness as well as a bending stiffness that acts transverse to the 
fiber contour. Previous studies have shown that this model displays an anomalous elastic regime where the 
stretching mode is fully frozen out and the elastic energy is completely dominated by the bending mode. 
We demonstrate by simulations and scaling arguments that, in contrast to the bending dominated elastic 
energy, the equally important elastic forces axe to a large extent stretching dominated. By characterizing 
these forces on microscopic, mesoscopic and macroscopic scales we find two mechanisms of how forces 
are transmitted in the network. While forces smaller than a threshold F c are effectively balanced by 
a homogeneous background medium, forces larger than F c are found to be heterogeneously distributed 
throughout the sample, giving rise to highly localized force-chains known from granular media. 
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1 Introduction 

It has been well known for more than a century that net- 
1 works of central force springs lose their rigidity when the 
number of connected neighbours is lower than a certain 
threshold value pQ . To guarantee the rigidity of these oth- 
' erwise "floppy" networks, additional contributions to the 
elastic energy, beyond central-force stretching, have to be 
introduced [2]- Here our focus is on a particular class of 
heterogeneous networks composed of crosslinked fibers, 
whose length If is much larger than the typical distance 
l s between two fiber-fiber intersections (see Fig.QJ. These 
systems have recently been suggested as model systems 
for studying the mechanical properties of paper sheets [3] 
or biological networks of semifiexible polymers 4,5,6 . As 
only two fibers may intersect at a given cross-link, the 
average number of neighboring cross-links is z < 4. This 
places them below the rigidity transition both in two and 
in three spatial dimensions. Several strategies have been 
used to elastically stabilize a central-force fiber network [7] . 
Here, we use an additional energy cost for fiber "bending" . 
The bending mode penalizes deformations transverse to 
the contour of the fiber while to linear order the distance 
between cross-links, i.e. the length of the fiber, remains 
unchanged. 

The two-dimensional fiber network we consider is de- 
fined by randomly placing N initially straight elastic fibers 
of length !j on a plane of area A — L? such that both po- 
sition and orientation are uniformly distributed. The fiber 
density is thus defined as p = Nlf/A. We consider the 
fiber-fiber intersections to be perfectly rigid, but freely 



rotatable cross-links that do not allow for relative sliding 
of the fibers. The elastic building blocks of the network 
are the fiber segments, which connect two neighbouring 
cross-links. A segment of length l s is modeled as a clas- 
sical beam with cross-section radius r and bending rigid- 
ity k [16]. Loaded along its axis (central force "stretch- 
ing") such a slender rod has a rather high stiffness, char- 
acterized by the spring constant fc||(^ s ) ~ n/l s r 2 , while 
it is much softer with respect to transverse deformations 
k±(ls) ~ k/1 3 s ("bending"). 

While this random fiber network is known to have a 
rigidity percolation transition at a density p c [HIKE!] i we 
have recently shown [TTl[T2j how the network's inherent 
fragility, induced by its low connectivity z, determines 
the properties even in the high-density regime far away 
from the percolation threshold, p 3> p c - In particular, it 
was possible to explain the anomalous scaling properties 
of the shear modulus as found in simulations [91ll3|. The 
unusually strong density dependence of the elastic shear 
modulus G ~ (p — pc) 6 67 [9] is found to be a consequence 
of the architecture of the network that features various 
different length-scales [TT | fT2 ) (see Fig. [TJ. On the meso- 
scopic scale the fiber length I f induces a highly non-affinc 
deformation field, where segment deformations <5 na follow 
the macroscopic strain 7 as 5 na ~ This is in stark 
contrast to an affine deformation field where deformations 
scale with the size of the object under consideration. For 
a segment of length l s the affine deformation therefore is 
Srmaff oc lh- Microscopically, a second length Z m ; n gov- 
erns the coupling of a fiber segment to its neighbours on 
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Fig. 1. fllustration of the local network structure and the rel- 
evant length-scales in the random fiber network (drawn in the 
deformed configuration): the fiber length If, the typical seg- 
ment length l a and the shortest deformed segment of length 
Imin- -Kcxt signifies the bending force that the crossing fiber ex- 
erts in the axial direction of the horizontal fiber, io and 
correspond to axial forces that are directly transmitted from 
one fiber to a neighboring fiber at the crosslink. This mecha- 
nism forms the basis for the establishment of force-chains; see 
main text. 

the crossing fiber [TT1I12) . Due to the fact that the bend- 
ing stiffness k± ~ l~ 3 of a segment is strongly increas- 
ing with shortening its length, it is found that segments 
with l s < ^ m j n rather deform their neighbours than being 
deformed, while longer segments l s > l min are not stiff 
enough to induce deformations in the surrounding. Thus 
^min plays the role of a rigidity scale, below which seg- 
ments are stiff enough to remain undcformed. As a con- 
sequence the elastic properties of the fiber as a whole are 
not governed by the average segment l s , but rather by the 
smallest loaded segment i m i n . 

In the observed scaling regime the elastic modulus 
does originate exclusively in the bending of the individual 
fibers, and thus reflects the stabilising effect of this soft 
deformation mode. In contrast, stretching deformations in 
this non-afnne bending dominated regime can be assumed 
to be frozen out completely. In our previous articles 11,12 
we have dealt with the properties of the bending mode, 
and its implications on the elastic energy. Here we 
focus on the stretching mode, ku , and its role in the oc- 
curence of elastic forces. As the non-afhne bending regime 
is relevant for slender fibers with a small cross-section 
radius, r — > 0, it is characterized by a scale separation, 
fcii/fcj_ ~ r~ 2 —¥ oo, which assures that no stretching de- 
formations, <5||, occur. Thus, the fibers effectively behave 
as if they were inextensible bars. Closer inspection of the 
limiting process reveals, however, that the stretching de- 
formations tend to zero as 5ii ~ 1/fcii El This makes the 

1 This may be seen by considering the following simplified 
energy function, W = k±(S — <5y) 2 + k\\5?i, which represents 
a system of two springs connected in series. Minimizing the 
energy for fixed overall deformation 5, one finds S« = kxS/(k\\ + 
k±), which shows the reqired scaling S\\ ~ in the limit 
fcn — » oo. 



contribution to the total energy W s ~ k^5 2 ~ k^ 1 negli- 
gible, as required, but also implies that finite stretching 
forces F\\ will occur, 

Fn ~ knSn — > const . (1) 

Indeed, these forces, acting axially along the backbone of 
the fibers, are absolutely necessary to satisfy force-balance 
at the intersection of two fibers. With two fibers intersect- 
ing at a finite angle it is intuitively clear that a change in 
transverse force in one fiber has to be balanced by an axial 
force in the second. 

For thicker fibers with a larger cross-section radius r, 
a second elastic regime occurs, where the bending instead 
of the stretching mode is frozen out. This regime for- 
mally corresponds to the limit kn/k± ~ r~ 2 — > 0. It is 
characterized by stretching deformations of mainly affine 
character [T3]. The elastic shear modulus (as well as the 
Young's modulus) have been shown to depend linearly on 
density [9,15,14 , G ~ p, which is in striking contrast to 
the strong susceptiblity to density variations found in the 
non-afnne bending regime. 

In the following we will present results of simulations 
that characterize in detail the occuring axial forces in both 
elastic regimes, the non- affine bending as well as the affine 
stretching regime. In the simulations we subject the net- 
work to a macroscopic deformation and determine the new 
equilibrium configuration by minimizing the elastic en- 
ergy. The minimization procedure is performed with the 
commercially available finite element solver MSC.MARC. 
Further details of the simulation procedure can be found 
in our previous publications ,6] . Starting with the average 
force-profile along the fibers we then proceed by giving 
the full probability distribution of forces. We show that 
its tails are heterogeneously distributed throughout the 
system, similar to force-chains in granular media. We find 
that most of the features can be understood with the help 
of the two basic length-scales, the filament length If and 
the rigidity scale Z m i n - 



2 Effective Medium Theory 

In this section we will characterize the configurationally 
averaged force profile along a fiber. In the spirit of an ef- 
fective medium theory, one can think of the fiber as being 
imbedded in an elastic matrix that, on a coarse-grained 
scale, acts continuously along the backbone. The associ- 
ated force /Coxt, which is imposed on the fiber, leads to a 
change in axial force according to [TB] 

dF \\ 

— — = -lC eKt , (2) 

OS 

where s is the arclength along the fiber backbone. 

A while back, Cox [T7] provided a second, constitutive 
equation that allows to solve for the force profile -Fj|(s). 
He assumed the medium to be characterized by an affine 
deformation field (W(s) ~ js. The external force /C ex t is 
assumed to be non- vanishing only when the actual fiber 
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deformation 8\i is different from this affine deformation 
field, 

/Cext(s) = k (S\\ - 5 af£ (s)) . (3) 

Eqs. ([2]) and ([3]) can easily be solved and result in 
a force profile that shows a plateau in the center of the 
fiber as well as boundary layers where the force decreases 
exponentially, Fn(s) = a — 6cosh[c(s — lf/2)], with a, b, 
and c appropriately chosen constants. 

Astrom et al. [H] have applied this model to the affine 
stretching regime and found the boundary layers to grow 
as the fiber cross-section radius is decreased. Fig. [2] shows 
results of our simulations for the force profile both in the 
affine stretching regime (blue squares) and the non-affinc 
bending regime (red circles). Apparently the Cox- model 
accounts very well for the force profile in the stretching 
regime, while it fails completely in the bending regime, 
where the simulation data clearly show that the force in- 
creases linearly from the boundary towards the center of 
the fiber. 
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Fig. 2. Variation of average axial force Fn along the backbone 
s — [0,lf] of the fibers (symmetrized around s = lf/2). The 
symbols are the results of our simulations. Towards the fiber 
ends the force relaxes exponentially in the affine stretching 
regime (full curve is a fit to the Cox model), and linearly in 
the non-affine bending regime. 

Cox ideas can be generalized to the non-affine bend- 
ing regime, where the elastic medium entirely consists of 
bending modes. In this regime the axial forces in the fiber 
arise due to the pulling and pushing of its crossing neigh- 
bors that try to transfer their high bending load in a kind 
of lever- arm effect (see Fig [J). As explained above, the de- 
formation field is non-affine and, instead of 6 a s, one has to 
use S na ~ jlf [123 , which is proportional to the fiber length 
//. Since we are interested in the limit where stretching de- 
formations vanish, Sn — > 0, the resulting exernal force is 
arclength independent, /C cx t = — kS na _ ~ — If and constant 
along the backbone. The axial force profile F«(s) is thus 
expected to be linearly increasing from the boundaries to- 
wards the center, in agreement with the simulations. 

Recently, Head et al. [13] have suggested growing bound- 
ary layers to play a key role for the cross-over from the 



affine stretching to the non-affine bending regime. Here, 
we have shown that these growing boundary layers are 
rather a consequence of a transition from an exponential 
to a linear force profile in the boundary layers. This follows 
naturally from the fact that non-affine deformations, <5 na , 
scale with the fiber length If and not with the segment 
length l s . As we have analyzed in detail in Ref. [12], this 
scaling property, which originates in the network architec- 
ture, can be understood within a "floppy mode" concept. 
Therefore, the growing boundary layers should be viewed 
as a consequence rather than the driving force of the affine 
to non-affine transition. 



3 Force distribution 

We now turn to a discussion of the full probability dis- 
tribution of axial forces, instead of just the average value 
as we have done in the previous section. In Fig. [3] we dis- 
play the distribution function P{Fii) for various densities 
p deep in the non-affine regime. Remarkably, the curves 
for different densities collapse on a single master curve by 
using the scaling ansatz 

P(F l{ )=F- l h(F {l /F c ), (4) 
with the force scale F c = Kp\{6p/ p c ) 5 , where Sp = p — p c . 
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Fig. 3. Probability distribution P{F^) (inset) and scaling 
function h in the non-affine bending regime with aspect ratio 
r/lf = 5 ■ 10" 6 for various densities plf = 20. ..100. The force 
scale used to obtain the data collapse is F c — npc(Sp/ p c ) 5 . 

Its appearance in Eq. ^ indicates that F c is the aver- 
age axial force. We now show that it is furthermore equiv- 
alent to the average bending force F c = (k±(l s ) ■ <5 na ) that 
is needed to impose the non-affine bending deformation 
<5 na on a segment of bending stiffness k± . 

The average (.) is performed over the segment-length 
distribution P(Z s )q restricted to those segments that are 

2 In the random network structure considered here, this dis- 
tribution is exponential, P(l s ) = pexp(— hp)- 
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longer than the rigidity scale Z m i n - Recall, that Z m ; n plays 
the role of the shortest deformed segment such that shorter 
segments do not contribute to the averaging. The equiva- 
lence of both expressions for the force scale becomes evi- 
dent by writing the average explicitly, 

/>OC 

F c = jf dl s P(l s )k ± (l s )S n& , (5) 

and inserting k± ~ l~ 3 . The integral is dominated by its 
lower limit, which leads to F c ~ n(pl f) 5 /I'j as required. 

Hence, the force scale F c has been identified as the av- 
erage force that induces the non-affine bending of a fiber. 
At the same time it occurs in the probability distribution 
of stretching forces as depicted in Fig. [31 This two-fold role 
reflects the interplay of bending and stretching in the ef- 
fective medium theory, where stretching forces in one fiber 
have to be balanced by bending forces in its neighbors. We 
can thus conclude that the effective medium picture is the 
appropriate description for forces up to the threshold F c , 
i.e. for the "typical" properties of the system. 

Interestingly, for forces smaller than F c the probabil- 
ity distribution displays an intermediate power-law tail 
h(x) ~ x~P, with an exponent (3 « 0.55, which as yet 
defies a simple explanation. For forces larger than F c , one 
may even speculate about the existence of a second power- 
law regime with a much higher exponent 0' « 5. While 
in this regime the distribution does not seem to decrease 
exponentially, the available range of forces is too small to 
reach any final conclusions as to the functional form. 

The axial forces in the affine stretching regime follow 
a completely different probability distribution, as can be 
seen from Fig. |4] The solid black line relates to the proba- 
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Fig. 4. Probability distribution P(-Fji) of axial forces in the 
affine stretching dominated regime for plf = 60 and aspect 
ratio r/lf — 0.05. The system is subject to uniaxial stretching, 
as indicated in the inset. Note the linear scale on the x-eods. 
The peak of the distribution is strongly correlated with the 
orientation, <f>, of the fiber relative to the principal stretching 
direction. 

bility distribution of all segments, irrespective of their ori- 
entation <fr, with respect to the imposed strain field (in this 



case: uniaxial extension along the j/-direction). The broken 
lines only include segments with orientations in a par- 
ticular interval, as denoted in the figure legend. The tails 
of the distribution are exponential (as has previously been 
observed in Ref. [14 ), while the peak force strongly de- 
pends on the orientation <f> of the segments. The position of 
the peak follows naturally from the assumption that seg- 
ments undergo affine deformations with S a g = jl s cos 2 4>. 
Since the resulting affine forces F a s((j)) = fcii(/ 8 )5aff ~ 
cos 2 4> do not depend on the length of the segments, one 
expects only a single, orientation-dependent force, that is 
a delta-function distribution P(F,<j>) ~ 5(F — F a f£ (</>)). 
The broadening of the distribution relative to the singu- 
lar delta-function can be rationalized with the fact that 
the affine strain field can only fulfill force equilibrium if 
the fibers are infinitely long. For finite fibers the force has 
to drop to zero at the ends (as discussed above), leading 
to additional (non-affine) deformations and therefore to a 
broadening of the peak. 



4 Force-chains 

In the preceding two sections we have characterized the oc- 
curing axial forces on the scale If of a single fiber as well as 
on the smaller scale of an indivdual fiber segment l s . We 
now proceed to discuss the properties of the forces on a 
larger, mesoscopic scale. To this end we probe the Green's 
function of the system and impose a localized perturba- 
tion in the center of the network. Fixing the boundaries, 
we displace a single crosslink and monitor the resulting 
response of the network. Similar studies have been con- 
ducted in Ref. [18], where ensemble averaging is used to 
discuss the applicability of homogeneous elasticity theory. 
Here, we focus on the individual network realization to 
better characterize the spatial organization of the forces. 

In Figs. [5] and [6] we show pictures of a network, where 
the grey-scales of the segments are chosen according to 
the magnitude of the forces they carry. The higher the 
force the darker the segment. While the quenched ran- 
dom structure is the same in both plots, the fiber aspect 
ratio a — r/lf is chosen such that the network lies deep 
in the non-affine bending (a = 10 -5 ) or the affine stretch- 
ing regime (a = 10 _1 ), respectively. As is clearly visible 
in Fig. [5] the non-affine response is characterized by well 
defined paths of high forces that extend from the center, 
where the force is applied, to the boundaries. These forces 
are transmitted from fiber to fiber along a zig-zag course. 
Upon comparison with the distribution function for axial 
forces (like those shown in Fig. [3]) , we find that the force- 
chains are constituted by forces with magnitude larger 
than F c (i.e. above the intermediate power-law regime). 
In contrast, in the affine regime a rather homogeneous 
distribution of forces is observed (Fig. [6J. 

The observation that the highest occuring forces are 
connected in chains, suggests a mechanism that allows di- 
rect transmission of an axial load from one fiber to the 
next. This has to be contrasted with the results of Sec.[2J 
where we have developed an effective medium picture that 
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Fig. 5. Picture of a network in the non-affine bending domi- 
nated regime. The grey-scale is chosen according to the axial 
force in the segment. Force-chains are clearly visible, and follow 
a zig-zag course from the center towards the lower boundary. 



Fig. 6. Picture of a network in the affine stretching dominated 
regime. The grey-scale is chosen according to the axial force in 
the segment. No force-chains are present. 



constitutes an indirect load transfer, in which axial forces 
in one fiber are balanced by transverse (bending) forces in 
the neighboring fiber. The difference between both mecha- 
nisms becomes evident by considering the load transfer at 
a fiber-fiber intersection that coincides with the end-point 
of one of the two fibers; see Fig. [7] for an illustration. An 
axial load Fu 1 ^ that is coupled into fiber (1) at its left 
end has to be taken up by the crossing fiber (2) at the 
distal end of fiber (1). Assume that this applied force in 




Fig. 7. Illustration of the two mechanisms of load transfer for 
a fiber segment (1) ending at a crossing fiber (2). The displace- 
ment &i of the horizontal fiber leads to deformations in both 
fibers. The indirect load transfer, b), couples axial forces into 
bending deformations, while the direct transfer, c), transmits 
the axial force directly into axial forces in the neighboring fiber. 

the horizontal fiber, iC , is accompanied with an axial 
displacement tiz\, which translates the fiber along its axis. 

In the indirect load transfer the axial force F^ 1%> is bal- 

(1 2) 

anced by bending forces Fj_ ' oriented perpendicular to 
the contour of both fibers. The force balance reads 

F^+F^=-F^ . (6) 

The angle between both fibers being 9, the associated 
bending displacements are (Fig. \7)p>) 

yx = —5zcot6, j/2 = <fe/sin#. (7) 

This mechanism forms the basis of the effective medium 
theory applied in Sec. [2] to explain the average force pro- 
file along the fibers. It has also been shown to allow a 
direct calculation of the macroscopic elastic moduli of the 
network [Tfl[T2]. 

We now proceed to discuss the second mechanism. It 
can immmediately be seen from the ^-dependence that the 
bending displacements j/1/2 of Eq.([7|) can become increas- 
ingly large, if the fibers intersect at small enough angles, 
8^0. Then, the network responds in a different way. 
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As can be inferred from Fig. [Th the large bending contri- 
butions may be avoided by additionally moving the sec- 
ondary fiber along its own axial direction, &z% = 8z\ cos 9. 
This leaves the primary fiber straight while the bending 
displacement in the secondary fiber becomes 1/2 — 5z\ sin 9. 
The force balance must therefore be written as 

F^ = -(Fp+F^). (8) 

The smaller the angle 9 the larger is the contribution of 

(2) 

Fn , which is just the forward "scattering" offerees, seen 
in Fig. [5] 

The existence of force-chains is intimately connected 
to the presence of long fibers, plf 3> 1. Since the coordi- 
nation can be written as z — z c (l — 0(plf)~ 1 ) [12], the 
system will become completely rigid by increasing the fiber 
length, If — > 00. There, the force scale F c diverges and ax- 
ial forces imposed on a fiber will propagate along the fiber 
completely uncorrelated with forces on its neighbours. As 
the number of intersections n c \ ~ pit per fiber is pro- 
portional to the fiber length, the probability distribution 
for the angle between two intersecting fibers is sampled 
to an increasing degree. In particular, the smallest angle 
that will occur in a finite sample of n c ; intersections is 
calculated from 




as m i n ~ If , where we have used P{9) ~ 9 for small 
values of 9. With increasing fiber length ever smaller an- 
gles occur. As explained above, the presence of small an- 
gles necessarily lead to forward scattering of axial forces, 
and thus to the emergence of force-chains. The presence 
of force-chains is therefore a consequence of the special 
geometry of the network that allows the fibers to intersect 
at arbitrarily small angles. To support this view, we have 
conducted additional simulations in which the localized 
perturbation is applied in varying directions. As a result, 
the structure of the force-chains remain unchanged while 
their amplitudes are modulated. 



5 Conclusion and Outlook 

In conclusion, we have characterized in detail the prop- 
erties of forces occuring in two-dimensional random fiber 
networks. We have shown that the previously identified 
rigidity scale Z m i n , in addition to the structural scale of 
the fiber length If, governs the occurence of stretching 
forces in an elastic regime, where the energy derives from 
the bending mode only (kti/kx — * 00). The probability 
distribution of forces shows scaling behaviour with the 
force scale F c = Kpl(5p/p c ) 5 that can be identified with 
the average force needed to deform a fiber segment by the 
non-affine deformation <5 na ~ 7Z/. 

Two types of force transmission have been identified. 
Forces up to the scale F c are transmitted from one fiber 
to the next by an indirect mechanism, where stretching 



forces are balanced by bending forces and vice versa. This 
is best illustrated by the action of a lever-arm that tries to 
transmit its bending load to the neighboring fiber, which 
subsequently starts to stretch (see /C GX t in FigQ]). This 
mechanism of force transmission can be used to under- 
stand the average force profile along a fiber, and also forms 
the basis for the calculation of the scaling properties of 
the elastic modulus [I1.12J. In a second direct mechanism 
axial forces are also transmitted directly to their neigh- 
bors, giving rise to highly localized force-chains that are 
heterogeneously distributed throughout the sample. This 
mechanism is only active for forces larger than F c . In con- 
trast to the indirect mechanism, which probes the center 
of the force distribution and therefore typical properties of 
the network, the direct mechanism reflects the properties 
of the extremal values of the distribution. 

The observation of force-chains suggests an analogy 
with granular media |19j . The scale separation ku 3> k± 
implies that the fibers behave as if they were inextensi- 
ble. The network of inextensible segments may therefore 
loosely be viewed as the contact network of a system of 
rigid grains. Due to the low coordination number, the fiber 
network has to be stabilised by the action of the bend- 
ing mode directed perpendicularly to the fiber axis. This 
bears some similarity to a friction force in granular sys- 
tems, which is directed tangentially to the grain surfaces. 
Indeed, it has been argued that stable granular systems 
with a coordination as low as found here may only be 
achieved if frictional forces are taken into account (20] . 
The occurence of force chains in our "frictional" system 
is, however, due to the vicinity of the isostatic point with 
regard to the frictionless system (z c = 4). Since the co- 
ordination may be written as z = 4(1 — O(plf)^ 1 ) the 
isostatic point is reached by increasing the fiber length 
If -► 00. 

We gratefully acknowldege fruitful discussions with Klaus 
Kroy. Financial support of the German Science Founda- 
tion (SFB486) and of the German Excellence Initiative 
via the program "Nanosystems Initiative Munich (NIM)" 
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